Analysis of scRNA-seq data from Pancreatic Islets
Abhijeet R. Patil
4/26/2022
Section 1- Preliminaries
Load libraries
To begin with the analysis, we will first load the required packages. The downstream analysis of the scRNA-seq is built upon Seurat and sctransform packages. The scDblFinder package is used for identifying and removing doublets in the dataset. The single-cell annotation is carried out using the marker-based scSorter method.
library("Seurat")
library("ggplot2")
library("sctransform")
library("BiocParallel")
library("scDblFinder")
library("SingleCellExperiment")
library("cowplot")
theme_set(theme_cowplot())
library("RColorBrewer")
library("viridis")
library("scSorter")
library("knitr")Load data
First, we will load the filtered count matrix from cellranger. This matrix contains the scRNA-seq counts from all the HPAP samples.
panc_data <- Read10X(data.dir = "./data/HPAP_samples_aggr_18Apr2022/outs/count/filtered_feature_bc_matrix/")Dimension of raw merged counts matrix
The dimension of the raw counts matrix after merging all HPAP scRNA-seq samples with cellranger is shown below. There are total of 36601 genes and 269417 cells.
# panc_data <- Read10X(data.dir = "PATH_TO_FEATURE_MATRIX")
dim(panc_data)## [1] 36601 269417Section 2- Preprocessing the data
Create Seurat object
The Seurat object initialization step above only considered cells that expressed at least 200 genes. The min.cells- parameter will include features detected in at least three cells. The min.features - Include cells where at least three fifty features are detected.
With the above parameters the count matrix will also be subsetted. Optionally, to revive the excluded features, create a new object with a lower threshold.
panc <- CreateSeuratObject(counts = panc_data, min.cells = 3, min.features = 200,
project = "hpap@67")
dim(panc)## [1] 32066 258379In the Seurat object, 32066 genes and 258379 cells were present after passing the thresholds.
Add metadata
Here, we will add four new columns to the Seurat metadata. The first column sample_id contains raw sample names from cellranger, hpap_id refers to the donor identifier, the disease_id corresponds to the numeric count in each disease type, and finally the disease_state shows the general state of disease (CTL, AAB, and T1D).
## Add metadata information
sample <- read.csv(file.path("./metadata", "HPAP_samples_aggr_18Apr2022_cellranger.csv"), stringsAsFactors=F)## Add sample_id
CellsMeta["sample_id"] <- df$sample_id
CellsMetaTrim <- subset(CellsMeta, select = c("sample_id"))
panc <- AddMetaData(panc, CellsMetaTrim)
## Add hpap_id
CellsMeta["hpap_id"] <- df$hpap_id
CellsMetaTrim <- subset(CellsMeta, select = c("hpap_id"))
panc <- AddMetaData(panc, CellsMetaTrim)
## Add disease_id
CellsMeta["disease_id"] <- df$disease_id
CellsMetaTrim <- subset(CellsMeta, select = c("disease_id"))
panc <- AddMetaData(panc, CellsMetaTrim)
## Add disease_state
CellsMeta["disease_state"] <- df$disease_state
CellsMetaTrim <- subset(CellsMeta, select = c("disease_state"))
panc <- AddMetaData(panc, CellsMetaTrim)Table showing metadata information
head(panc@meta.data) %>%
kable(format = "html", col.names = colnames(head(panc@meta.data))) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "100px")| orig.ident | nCount_RNA | nFeature_RNA | sample_id | hpap_id | disease_id | disease_state | |
|---|---|---|---|---|---|---|---|
| AAACCTGAGCGTTGCC-1 | hpap@67 | 16875 | 4080 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGAGGACAGAA-1 | hpap@67 | 6497 | 892 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGAGTCGTACT-1 | hpap@67 | 1169 | 253 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGCAACTGCTA-1 | hpap@67 | 13689 | 3553 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGCAATGTAAG-1 | hpap@67 | 16375 | 3583 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
| AAACCTGCACCAGGTC-1 | hpap@67 | 16431 | 3438 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D |
## Calculate MT %
panc[["percent.mt"]] <- PercentageFeatureSet(panc, pattern = "^MT-")The below plot shows the number of features, counts, and mitochondrial reads before filtering.
VlnPlot(panc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0.001, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
VlnPlot(panc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
Correlation between nCount_RNA - nFeature_RNA and nCount_RNA - percent.mt before filtering.
# FeatureScatter used to visualize feature-feature relationships.
plot1 <- FeatureScatter(panc, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "disease_state") & theme(legend.position = 'none')
plot2 <- FeatureScatter(panc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "disease_state")
CombinePlots(plots = list(plot1, plot2))
Remove doublets
The scDblFinder function is employed to remove the potential doublets from each HPAP donor separately. The parameter sample_id refers to individual donor. The function uses an object from SingleCellExperiment class, therefore the SingleCellExperiment package is used for conversion.
## Convert object into singlecellexperiment
panc.sce <- as.SingleCellExperiment(panc)
panc.sce <- scDblFinder(panc.sce, samples="sample_id", clusters=FALSE, BPPARAM=MulticoreParam(10))
## Convert sce object back to seurat
panc_seurat <- as.Seurat(panc.sce, counts = "counts", data = "logcounts")The total number of doublets are shown below
table(panc.sce$scDblFinder.class)##
## singlet doublet
## 234285 24094The doublet information is stored in the metadata. The updated metadata can be found below
Updated metadata table with doublets information
head(panc@meta.data) %>%
kable(format = "html", col.names = colnames(head(panc@meta.data))) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "100px")| orig.ident | nCount_RNA | nFeature_RNA | sample_id | hpap_id | disease_id | disease_state | percent.mt | |
|---|---|---|---|---|---|---|---|---|
| AAACCTGAGCGTTGCC-1 | hpap@67 | 16875 | 4080 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 7.330370 |
| AAACCTGAGGACAGAA-1 | hpap@67 | 6497 | 892 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 1.062029 |
| AAACCTGAGTCGTACT-1 | hpap@67 | 1169 | 253 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 4.362703 |
| AAACCTGCAACTGCTA-1 | hpap@67 | 13689 | 3553 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 3.382278 |
| AAACCTGCAATGTAAG-1 | hpap@67 | 16375 | 3583 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 3.877863 |
| AAACCTGCACCAGGTC-1 | hpap@67 | 16431 | 3438 | HPAP021_75162_T1DIslets | HPAP021 | T1D_ 1 | T1D | 1.545858 |
Pre-processing the data based on information from above steps
The nFeature_RNA and nCount_RNA represents the number of genes detected in each cell and the total number of molecules (UMIs) detected within a cell respectively. While the low nFeature_RNA in a cell mean that it may be dead/dying or it may represent an empty droplet, the high nCount_RNA and nFeature_RNA denotes that the cell may be a doublet or multiplet. This filtering along with mitochondrial reads is crucial pre-processing step, because, removing such outliers from these groups might also remove some of the doublets or dead/empty droplets.
Remove all the potential doublets found
The cell being either doublet or singlet is stored in column scDlFinder.class. The potential doublets are removed in this step.
# After doublet removal
panc <- subset(panc_seurat, subset = scDblFinder.class == "singlet")The data after filtering based on nFeature and percent.mt is shown below
panc_new <- subset(panc, subset = nFeature_RNA > 200 & nFeature_RNA < 9000
& percent.mt < 25
& nCount_RNA < 120000
)
dim(panc_new)## [1] 32066 222679Update Metadata with MT and Rb genes information
After removing doublets, now we can proceed with further pre-processing steps. For ex, Mitochondrial genes are useful indicators of cell state.
Add the mitochondrial proteins information to the metadata after fitering
Adding the Mitochondrial genes
panc[["percent.mt"]] <- PercentageFeatureSet(panc, pattern = "^MT-")Display the metadata information after filtering
The below plot shows the number of features, counts, and the percentage of mitochondrial reads after filtering.
VlnPlot(panc_new, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0.001, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
VlnPlot(panc_new, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = "disease_state",
pt.size = 0, combine = T) & theme(legend.position = 'none',
axis.title.x = element_blank())
Correlation between nCount_RNA - nFeature_RNA and nCount_RNA - percent.mt after filtering.
The FeatureScatter is typically used to visualize feature-feature relationships.
plot1 <- FeatureScatter(panc_new, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "disease_state") & theme(legend.position = 'none')
plot2 <- FeatureScatter(panc_new, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "disease_state")
CombinePlots(plots = list(plot1, plot2))
Section 3- Normalizing and scaling the data
SCT: Single Cell Transform
The seurat’s SCTranform (SCT) function is used to normalize the counts that measures the differences in sequencing depth per cell for each sample. The SCT method is built on the concept regularized negative binomial model to perform normalization and variance stabilization of single-cell RNA-seq data. It removes the variation due to sequencing depth (nUMIs). The vars.to.regress parameter allows us to regress out the variation from other sources such as percentage of mitochondrial. The output of SCT model is the normalized expression levels for all the transcripts. Lastly, the variable.features.n parameter is used to select the variable features and is set as 3000.
panc_new_sct <- SCTransform(panc_new, method = "glmGamPoi", variable.features.n = 3000, vars.to.regress = "percent.mt", verbose = FALSE)The seurat object now consists of two assays RNA and SCT, with the default being set as SCT.
DefaultAssay(panc_new_sct) ## [1] "SCT"Perform dimensionality reduction by PCA and UMAP embedding
These are now standard steps in the Seurat workflow for visualization and clustering
panc_new_sct <- RunPCA(panc_new_sct, verbose = TRUE)Plot showing heatmaps of PCs- 1 to 12 from the PCA analysis
DimHeatmap(panc_new_sct, dims = 1:12, cells = 500, balanced = TRUE)
Plot showing heatmaps of PCs- 13 to 24 from the PCA analysis
DimHeatmap(panc_new_sct, dims = 13:24, cells = 500, balanced = TRUE)
Plot showing heatmaps of PCs- 25 to 36 from the PCA analysis
DimHeatmap(panc_new_sct, dims = 25:36, cells = 500, balanced = TRUE)
Plot showing the variance explained by PCs from the PCA analysis
ElbowPlot(panc_new_sct, ndims = 40)
ElbowPlot(panc_new_sct)
Run UMAP
The Uniform Manifold Approximation and Projection (UMAP) dimensional reduction technique is run on the first 20 principal components.
panc_new_sct <- RunUMAP(panc_new_sct, dims = 1:20, verbose = FALSE)Find Neighbors
Here we compute the k nearest neighbors in the dataset. The nearest neighbors are calculated on the reduced dimensions from the PCA. We use first 20 principal components.
panc_new_sct <- FindNeighbors(panc_new_sct, dims = 1:20, verbose = FALSE)Find Clusters
Once we have calculated the k-nearest neighbors in previous step, we will find the clusters of cells by the shared nearest neihbor (SNN) clustering algorithm. Here, we have set the resolution as 0.8. The higher resolution usually gives more number of clusters.
panc_new_sct <- FindClusters(panc_new_sct, resolution = 0.8, method = "igraph", verbose = FALSE)UMAP plotting
UMAP plot showing the clusters based on SCT data and snn optimization at resolution 0.8
DimPlot(panc_new_sct, group.by='SCT_snn_res.0.8', reduction='umap') +
ggtitle('SCT_snn_res.0.8')
UMAP plot showing the distribution of cells across different disease types.
DimPlot(panc_new_sct, group.by='disease_state', reduction='umap') +
ggtitle('Disease type')
UMAP plot showing the cell distribution with respect to individual donors
DimPlot(panc_new_sct, group.by='hpap_id', reduction='umap') +
ggtitle('Sample type') 
Section 4- Cell Annotation
scSorter: single cell Sorter
Loading marker genes data
anno <- read.csv("./annot/anno.csv")
anno[,1] <- NULLDisplay marker genes across celltypes
We used marker-based scSorter cell annotation method for assigning cells to known cell type based on the marker genes. The markers for each cell type were added based on comprehensive literature search.
anno %>%
kable(format = "html", col.names = colnames(anno)) %>%
kable_styling() %>%
kableExtra::scroll_box(width = "100%", height = "100px")| Type | Marker |
|---|---|
| Acinar | PRSS1 |
| Acinar | REG1A |
| Acinar | CPA1 |
| Acinar | CPA2 |
| Alpha | GCG |
| Alpha | GC |
| Alpha | TTR |
| Beta | INS |
| Beta | IGF2 |
| Beta | IAPP |
| Ductal | KRT19 |
| Ductal | CFTR |
| Ductal | SFRP5 |
| Ductal | MMP7 |
| Stellates_Mesenchymal | COL1A1 |
| Stellates_Mesenchymal | PDGFRB |
| Stellates_Mesenchymal | RGS10 |
| Stellates_Mesenchymal | THY1 |
| Endothelial | VWF |
| Endothelial | CD93 |
| Immune | NCF2 |
| Immune | PTPRC |
| PP_Gamma | PPY |
| PP_Gamma | CARTPT |
| PP_Gamma | PCDH10 |
| PP_Gamma | PLAC8 |
| Delta | SST |
| Delta | LEPR |
| Delta | PRG4 |
| Epsilon | GHRL |
The scSorter method is based on a semi-supervised learning algorithm. It expects the pre-processed data as input. Therefore, we used the previously normalized, scaled and transformed expression matrix generated by SCT method as input. The top 3000 highly variable genes selected by SCT method was also given as input for scSorter method.
expr <- panc_new_sct@assays$SCT@data
topgenes <- head(panc_new_sct@assays$SCT@var.features, 3000)
topgene_filter = rowSums(expr[topgenes, ]!=0) > ncol(expr)*.1
topgenes = topgenes[topgene_filter]
## At last, we subset the preprocessed expression data and run scSorter.
picked_genes = unique(c(anno$Marker, topgenes))
expr = expr[rownames(expr) %in% picked_genes, ]fit the scSorter method
Based on the expression matrix and the marker genes the scSorter method is fit and the cell type assignment results predicted are stored in the Pred_Type vector.
rts <- scSorter(expr, anno)
panc_new_sct$cell_type <- rts$Pred_Type
table(rts$Pred_Type)##
## Acinar Alpha Beta
## 59496 65075 45009
## Delta Ductal Endothelial
## 5664 20300 5864
## Epsilon Immune PP_Gamma
## 196 1220 2821
## Stellates_Mesenchymal Unknown
## 9583 7451Add celltype information into the metadata
The celltype information is added to the metadata in the original Seurat object
cell_type.info <- data.frame(cell_type = panc_new_sct$cell_type, row.names= colnames(panc_new_sct_umap))
panc_new_sct_umap <- AddMetaData(object = panc_new_sct_umap, metadata = cell_type.info)
## Convert cell type to factor
panc_new_sct_umap$cell_type <- as.factor(panc_new_sct_umap$cell_type)
Idents(panc_new_sct_umap) <- "cell_type"
panc_new_sct_umap$ident <- NULLCheck expression of markers
The expression of markers across all the cell types are shown in different UMAPs
FeaturePlot(panc_new_sct_umap, features = c("PRSS1", "REG1A", "CPA1", "CPA2", "GCG",
"GC", "TTR", "INS", "IAPP", "GHRL",
"COL1A1", "PDGFRB", "RGS10", "THY1", "VWF",
"CD93", "PPY", "SST", "NCF2", "PTPRC"), pt.size = 0.2, ncol = 5) + theme(legend.position = 'none',
axis.title.x = element_blank(), axis.title.y = element_blank())
The UMAP showing the cell type annotations from the scSorter method
DimPlot(panc_new_sct_umap, group.by = "cell_type", reduction = "umap", label = FALSE, cols= c("#e30800", "#f56505", "#dec400", "#006630", "#0223c7","#5b02c7", "#00b0e6", "#c40080", "#02f00a", "#7d3301", "#000000"))
Below we show the UMAPs for the top 5 cell types
UMAP showing only Acinar cells
DimPlot(panc_new_sct_umap, group.by = "cell_type", reduction = "umap", label = FALSE, cols= c("#fac720", "#000000", "#000000", "#000000", "#000000","#000000", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000"))
UMAP showing only Alpha cells
UMAP showing only Beta cells
UMAP showing only Ductal cells
UMAP showing only Immune cells
